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ABSTRACT 


Energy storage components improve the energy efficiency of systems by reducing the mismatch 
between supply and demand. For this purpose, phase-change materials are particularly attractive since 
they provide a high-energy storage density at a constant temperature which corresponds to the phase 
transition temperature of the material. Nevertheless, the incorporation of phase-change materials 
(PCMs) in a particular application calls for an analysis that will enable the researcher to optimize 
performances of systems. Due to the non-linear nature of the problem, numerical analysis is generally 
required to obtain appropriate solutions for the thermal behavior of systems. Therefore, a large amount 
of research has been carried out on PCMs behavior predictions. The review will present models based on 
the first law and on the second law of thermodynamics. It shows selected results for several 
configurations, from numerous authors so as to enable one to start his/her research with an exhaustive 
overview of the subject. This overview stresses the need to match experimental investigations with 
recent numerical analyses since in recent years, models mostly rely on other models in their validation 
stages. 

© 2010 Published by Elsevier Ltd. 
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1. Introduction 


The ever increasing level of greenhouse gas emissions 
combined with the overall rise in fuel prices (although fluctuations 
occur) are the main reasons behind efforts devoted to improve the 
use of various sources of energy. Economists, scientists, and 
engineers throughout the world are in search of (1) strategies to 
reduce the demand; (2) methods to ensure the security of the 
supplies; (3) technologies to increase the energy efficiency of 
power systems; (4) new and renewable sources of energy to 
replace the limited and harmful fossil fuels. 

One of the options to improve energy efficiency is to develop 
energy storage devices and systems in order to reduce the 
mismatch between supply and demand. Such devices and systems 
also improve the performance and reliability by reducing peak 
loads and allowing systems to work within an optimal range. 
Thus, they play a preponderant role in conserving energy. The 
different forms of energy that can be stored are mechanical, 
electrical, and thermal. Here, mechanical (gravitational, com- 
pressed air, flywheels) and electrical (batteries) storages are not 
considered while thermal energy storage is discussed in the 
context of latent heat (sensible heat and thermochemical heat are 
not considered). 

Latent heat storage is based on the capture/release of energy 
when a material undergoes a phase change from solid to liquid, 
liquid to gas, or vice versa. Latent heat storage is particularly 
attractive since it provides a high-energy storage density and has 
the capacity to store energy at a constant temperature - or over a 
limited range of temperature variation - which is the temperature 
that corresponds to the phase transition temperature of the 
material. For instance, it takes 80 times as much energy to melt a 
given mass of water (ice) than to raise the same amount of water by 
1 °C. Table 1 provides a typical comparison between properties of 
different thermal storage materials used at room temperature. For 
the interested reader, excellent global reviews that pertain to 
phase-change materials and their various applications were 
proposed by Farid et al. [1], Sharma et al. [2], Zhang et al. [3], 
Regin et al. [4], Tyagi and Buddhi [5], Mondal [6], Sethi and Sharma 
[7], and especially the recent one by Verma et al. [8]. 

Nevertheless, the incorporation of phase-change materials 
(PCMs) in a particular application calls for an analysis that will 
enable the researcher to determine whether or not PCMs will 
improve performances sufficiently to justify extra costs for 
additional systems and/or controls needed. Mathematical model- 
ing of latent heat energy storage materials and/or systems is 
needed for optimal design and material selection. Therefore, a large 
amount of research has been carried out on PCMs behavior 
predictions whether they are considered separately or within 
specific systems. 


Table 1 shows the main differences, on average basis, between 
sensible heat storage materials and latent heat storage classes of 
materials. The table indicates that generally, PCMs require much 
less mass or volume to store the same amount of energy at a more 
or less constant temperature. 

This paper is building on previous reviews [1-8] to update the 
available references that pertain to mathematical modeling and 
simulation of thermal energy storage with phase-change materi- 
als. First, it presents the fundamental mathematical description of 
the phenomenon, the Stephan problem. Then, it provides basic 
mathematical descriptions used as basis for numerical modeling 
using either first or second law approaches and fixed or adaptative 
meshes. The next section, considered by the authors as the major 
contribution, is a model collection of most recent works published 
on the subject. This survey is organized according to the problem 
geometry (Cartesian, spherical, and cylindrical) and specific 
configurations or applications (packed beds, finned surfaces, 
porous and fibrous materials, slurries). A synthesis is provided 
at the end and several recommendations are formulated. 


2. The Stephan problem 


Phase transition of a material is described by a particular kind of 
boundary value problems for partial differential equations, where 
phase boundary can move with time. This question has been first 
studied by Clapeyron and Lamé in 1831 when analyzing the 
formation of the Earth’s crust when it cooled. In that case, the 
problem was simplified from a spherical geometry to a one- 
dimensional semi-infinite slab [9]. This solution was found 
independently by Franz Neumann, who introduced it in his 
lectures notes of 1835-1840 [10]. Nevertheless, this type of 
problems is named after Jozef Stefan, the Slovene physicist who 
introduced the general class of such problems in 1889 [11] in 
relation to problems of ice formation. Existence of a solution was 
proved by Evans in 1951 [12], while the uniqueness was proved by 
Douglas in 1957 [13]. 

Very few analytical solutions are available in closed form. They 
are mainly for the one-dimensional cases of an infinite or semi- 
infinite region with simple initial and boundary conditions and 
constant thermal properties. Under these conditions, these exact 
solutions usually take the form of functions of the single variable x/ 
t'/? and are known as similarity solutions [14,15]. A collection of 
similarity solutions and references is to be found in [16,17]. 


3. Numerical solution 
The problem of predicting the behavior of phase-change 


systems is difficult due to its inherent non-linear nature at moving 
interfaces, for which displacement rate is controlled by the latent 


Table 1 
Common heat storage materials. 
Property Materials 
Rock Water Organic PCM Inorganic PCM 
Density [kg/m?] 2240 1000 800 1600 
Specific heat [kJ/kg K] 1.0 4.2 2.0 2.0 
Latent heat [kJ/kg] - 190 230 
Latent heat [MJ/m?] = = 152 368 
Storage mass for 10°J, avg [kg] 67,000 16,000 5300 4350 
(At=15K) (At=15K) 
Storage volume for 10°J, avg [m°] 30 16 6.6 27 
(At=15K) (At=15K) 
Relative storage mass 15 4 1.25 1.0 
(At=15K) (At=15K) 
Relative storage volume 11 6 2.5 1.0 
(At=15K) (At=15K) 
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heat lost or absorbed at the boundary. The following equation, 
known as the Stephan condition, describes this process: 


(E = (EB) E) o 
where à is the latent heat of fusion, p is the density (it is not 
specified if it is solid or liquid), s(t) is the surface position, k is the 
thermal conductivity, t is time, and T is the temperature. Indexes s 
and | refers to solid and liquid phases. 

In this situation, the position and the velocity of boundaries are 
not known a priori. In addition, since two phases possess different 


physical properties, this could create the in numerical model non- 
physical discontinuities which need to be addressed. 


3.1. Fixed grid 


By introducing an enthalpy method, the phase-change problem 
becomes much simpler since the governing equation is the same 
for the two phases; interface conditions are automatically achieved 
and create a mushy zone between the two phases. This zone avoids 
sharp discontinuities that may create some numerical instabilities. 
In consequence, the thickness and the quality of the discretization 
of this mushy zone are critical to the model performance. The 
enthalpy method can deal with both mushy and isothermal phase- 
change problems but the temperature at a typical grid point may 
oscillate with time [18]. This method has been successfully applied 
to various phase-change problems [19-22]. Hunter in 1989 [23] 
and Amdjadi in 1990 [24] confirmed that the enthalpy method is 
the most suitable for typical applications under the restriction that 
there is no alteration to the numerical scheme at the interface. 

Enthalpy function h, defined as a function of temperature, is given 
by Voller [25]. For a phase-change process, energy conservation can 
be expressed in terms of total volumetric enthalpy and temperature 
for constant thermophysical properties, as follows (from [2]): 
oF L Viki(VT)), (2) 
where H is the total volumetric enthalpy, which is the sum of 
sensible and latent heats: 


H(T) = A(T) + py f(T)A (3) 
and where 

T 
h= f. PrCk dT. (4) 


In the case of isothermal phase change, the liquid fraction fis given 
by 


(0) T<Tm solid, 
f=4 ]0,1[ T=Tm mushy, (5) 
1 T>Tm liquid. 


Using Eqs. (3) and (4), one can write an alternate form for 
dimensional heat transfer in the PCM in two dimensions: 


oh a oh ð oh of 
BE ~ OX (o a Jy (« y) Pih ge (6) 
where a is the thermal diffusivity. 

Very early, two approaches, finite difference and finite-element 
techniques, were used to solve the phase-change problems 
numerically. Then, finite volume and control volume finite- 
element methods were also employed [26]. Since the introduction 
of finite-element and finite-volume methods in the 1970s, they are 
mostly preferred over finite-difference methods. 

Since temperature can oscillate with time in some circum- 
stances, some authors have relayed on a temperature based heat 


capacity method [27]. In these models, the enthalpy-based energy 
equation is converted into a non-linear equation with a single 
variable. This approach is called the temperature transforming 
model (TTM). This type of simulation has proved accurate, simple 
and efficient [27]. 

In TTM, general continuity and momentum equations for fluid 
problems are used. But unlike in the enthalpy formalism, the 
energy equation uses a temperature dependant source term. The 
main difficulty of this technique is to develop a method to keep 
velocity at zero in the solid phase. The simplest method is the 
switch-off method (SOM) or its smoothed version the ramped 
switch-off method (RSOM) [20,28,29]. These techniques directly set 
velocities in the solid phase to zero by setting coefficients of 
momentum and velocity-correction equations to values that 
effectively suppress any motion. Although these methods are 
commonly used in phase-change problems, Ma and Zhang [18] 
have shown that if used together with a TTM model for convection 
controlled solid-liquid phase-change problems, this will result ina 
serious inconsistency. This is because TTM uses a mushy region to 
guarantee that the temperature is continuous while in SOM, 
velocities are discontinuous at the solid—liquid interface. To avoid 
this discontinuity, a variant of this method uses a ramped switch- 
off of velocities by introducing a mushy zone between solid and 
liquid phases. 

Instead of modifying coefficients of momentum and velocity- 
correction equations, an alternate approach modifies the source 
term (STM) in the energy equation at the interface. As SOM, this 
method also presents a discontinuity at the interface and, 
therefore, suffers from the same problem. Therefore, a ramped 
version of this approach should also be used (here RSTM). It should 
be noted that Darcy STM developed for phase-change simulations 
in the context of enthalpy method is a type of ramped STM [20,30]. 
Finally, a third approach is to work with a variable viscosity of the 
medium (VVM) as proposed by Gartling [31]. In this method, 
viscosity is set to a very high value in the solid phase, which 
effectively stops all motions. 

Ma and Zhang [18] did a comparison between each of the TTM 
methods and the experimental results of Okada [32]. They 
concluded that ramped SOM and ramped STM performed better 
at a slightly higher computational cost. In addition, they pointed 
out that if the time step is too small, the convection is not well 
modeled and that the solutions of the model diverge. 


3.2. Adaptive mesh 


To keep a reasonable representation of the physical processes 
that occur at the melt front, model grid density must be high 
enough to smoothly cover the solid-liquid interface. However, 
this high density is not needed elsewhere in the numerical 
domain. Therefore, it is natural to adapt the grid density to the 
local physical conditions to improve the computational efficiency. 
There are two main approaches used to do this. One uses a local 
mesh refinement method, i.e. h-method [33-36]. In this case, the 
model starts with a uniform grid and, at each iteration, grid points 
are added or removed to match the required accuracy. This is the 
method used in most commercial codes. The main problem with 
this method consists in maintaining data structures since the 
topology and the number of grid points is changing between each 
time step. 

To avoid this problem, the r-method (r stand for relocation), 
also known as the moving mesh method, starts with a uniform mesh 
and then moves the mesh points, keeping the mesh topology and 
number of mesh points fixed as the solution evolves. Grid 
deformation is usually done by tracking a rapid variation of either 
the solution or one of its higher order derivatives. This method has 
been used with phase-change problems as in [37-40]. 
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Lacroix and Voller [41] performed a study comparing methods 
of simulation of a phase-change model in a rectangular cavity. They 
concluded that the fixed grid must be finer for a material with a 
unique melting temperature while the limiting factor in moving 
mesh is the need to use a coordinate generator at each time 
increment. Comparison between fixed and moving grid has also be 
done by Viswanath and Jaluria [42] and by Bertrand et al. [43]. In 
this last paper, it was found that front-tracking methods are better 
adapted to the problem than the fixed-grid procedures. The front- 
tracking methods would however fail to simulate situations where 
the transition from the liquid to the solid phase is not a 
macroscopic surface, and enthalpy methods are to be used in 
most solidification problems where a solid-liquid interfacial 
region is present between both phases. 


3.3. Fist law and second law models 


Some models are designed to measure the first law efficiency 
while others are for the second law. First law models have some 
shortcomings because they do not consider the effect of time 
duration through which the heat is stored or retrieved, the 
temperature at which the heat is supplied, and the temperature 
of the surroundings. Second law models address those issues which 
lead to optimal designs and operations. Nevertheless, second law 
models are intended to complement not replace the first law models. 

The first law efficiency of a thermal energy storage system is 
defined by 


-Qr 
=g; 


where Q, is the total thermal energy extracted from the storage 
material during the heat recovery process and Q, is the total heat 
stored by the phase-change material during the heat storage 
process. 

The second law efficiency of the thermal energy storage system 
is defined by the following equation: 


(7) 


pee (8) 


where @, is the availability recovered from the storage system during 
the heat recovery process and @, is the availability added to the 
system during the heat storage process. This shows that the first law 
efficiency gauges how energy is utilized, whereas the second law 
efficiency indicates how well the availability of energy is used. 

An important feature in the second law analysis of the thermal 
energy storage system is the calculation of the entropy generation 
number N,. This number is defined as the degraded exergy divided 
by the total input exergy to the cycle: 


Wa 


NS = WW, 


(9) 
where W, is the total exergy destroyed throughout the cycle, while 
W, and W, are the availability of the working fluid entering the 
system during the heat storage and removal processes, respective- 
ly. The second law efficiency can be related to the number of 
entropy generation units by the following equation: 


= 1—Nz. (10) 


In this situation, the system achieves its maximum second law 
efficiency when N, is zero. 


4. Model collection 


In this section, we review the scientific literature on numerical 
models of PCM used in thermal storage applications. This review is 


based on compilations by Verma et al. [8] and Regin et al. [4] 
combined with additional articles found trough a recent literature 
review. These papers have been first classified by their geometrical 
configuration for more convenience and then by applications. 


4.1. Rectangular geometry 


In a pioneer work, Shamsundar and Sparrow [44,45] applied the 
finite-difference method to the resolution of the enthalpy equation 
in the solidification of a flat plate. Shamsundar also used this 
method to the case of a square geometry [46-48]. 

Hamdam and Elwerr [49] presented a simple model where all 
sides of a rectangular enclosure were well insulated except one 
vertical side where heat was applied. In this model, the rate of 
melting depends essentially on the properties of the material, such 
as thermal diffusivity, viscosity, conductivity, latent heat of fusion 
and specific heat. The model follows a two-dimensional melting 
process of a solid PCM considering convection as the dominant 
mode within the melt region, except within the layer very close to 
the solid boundary in which only conduction is assumed to take 
place. Comparison between this model and results of Bernard et al. 
[50] shows a good agreement except that position and inclination 
of the melt front deviate as time progresses. 

Lacroix et al. [51-54] solved the problem of fusion in a 
rectangular cavity including the natural convection effects using a 
methodology similar to the moving mesh technique. Results of 
these simulations indicate that the melted fraction from close- 
contact melting at the bottom of the cavity is larger, by an order of 
magnitude, than that from the conduction dominated melting at 
the top. The melting process is essentially governed by the 
magnitude of the Stefan number. It is also strongly influenced by 
lateral dimensions of the cavity. The importance of the convection 
also implies that system efficiency is increased when the heat is 
concentrated on few spots instead of being spread over the whole 
surface. Later Brousseau and Lacroix [55,56] studied performances 
of a multi-layer PCM storage unit. Their model is based on the 
conservation equation of energy for the PCM and for the HTF (heat- 
transfer fluid). This model shows that cyclic melting and 
solidification yields complex isotherm distributions within the 
PCM with multiple phase fronts. 

Works by Costa et al. [57] studied numerically a two- 
dimensional rectangular area, using energy equations in solid 
and liquid phases, continuity, momentum, and Stefan’s equation at 
the boundary. They applied the enthalpy-porosity approach used 
by Voller et al. [30,250] to the governing equations. Three PCMs 
were analyzed: paraffin (n-octadecanol) and metals (gallium and 
tin). They compared their results with those obtained in the 
literature [245-249]. In the case of octadecanol, there is poor 
agreement with experimental results in the upper zone where melt 
liquid from the sides fills the upper empty cavity and accelerates 
melting in this area. Variation of viscosity with temperature, 
heating losses in the wall or different initial supercooling are 
proposed explanations for this discrepancy [57,58]. 

Since most PCMs have low thermal conductivities, then heat- 
transfer rates limit their applications. To improve their perfor- 
mances they are often used in a thin plate configuration like a heat 
exchanger [59-64]. Costa et al. [65] also studied a thermal storage 
system consisting of seven thin rectangular containers of PCM. The 
shape of the container is used to enhance the rate of heat transfer. 
The performance of this system has been analyzed using an 
enthalpy formulation and a fully implicit finite-difference method 
[25,66]. Vakilaltojjar and Saman [67] proposed to improve on this 
design by using PCM with different melting point (CaClz-6H20 and 
potassium fluoride tetrahydrate (KF-4H,0). 

Dolado et al. [68] developed a model for thermal energy storage 
unit for air conditioning application using thin PCM slab. In 
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parallel, they have also characterized the PCM in laboratory as well 
as experimentally validated their model. Comparing various 
numerical approaches, they showed that a simple 1D implicit 
finite differences model seems to be an appropriate model to 
simulate both, single PCM slabs and full storage systems. 
Comparing the advantages of various numerical approaches, they 
concluded that semi-analytical models give a fast first approach, 
successful to optimize the design, while 2D model should be used 
to study thicker plates and that the complete computational fluid 
dynamics model should be employed for more rigorous analyses. 

Zukowski [69] has analyzed the heat and mass transfer in a 
ventilation duct filled with encapsulated paraffin wax RII-56. This 
author proposed a new approach for approximating the specific heat 
of the PCM as a function of its temperature for the whole range of 
operating conditions, by using an interpolating cubic spline function 
to estimate specific heat of PCM at each temperature. His analysis 
showed that the average charge time (melting) of the tested units 
was over twice higher than the time of the discharge cycle (freezing). 

Silva et al. [70] have modeled an enclosure filled with wax in a 
rectangular geometry oriented vertically. Their results show that a 
simplified numerical model can be used to predict, with reasonable 
accuracy, the performance of this kind of heat exchange unit. 
Vynnycky and Kimura [71] produced an analytical and numerical 
study of coupled transient natural convection and solidification in 
a rectangular enclosure. Their work, based on a non-dimensional 
analysis, indicates that some asymptotic simplification is possible 
for materials commonly used in literature (water, gallium, lauric 
acid). This suggests that problems might be simplified, when 
Ra > 1 and St < 1, giving a conventional boundary value problem 
for the liquid phase and pointwise-in-space first-order ordinary 
differential equations for the evolution in time of the solidification 


front. The method was tested against full 2D finite-element-based 
transient numerical simulations of solidification. The asymptotic 
analysis decouple the fluid flow and heat-transfer problem in the 
liquid from the heat-transfer problem in the solid, and is able to 
describe the quantitative features of the numerical solutions very 
well for about 90% of the height of the enclosures. However, in the 
final 10%, the analytical solution is not uniformly valid for all time, 
and tends to overestimate the final thickness of the solid layer. 

Zivkovic and Fujii [72] have created a simple computational 
model for isothermal phase change. The model was based on an 
enthalpy formulation with equations constructed in such a form 
that the only unknown variable is the PCMs temperature. The 
theoretical model was verified with a test problem and an 
experiment performed by authors. Experimental and computa- 
tional data demonstrated that heat conduction within the PCM in 
the direction of the fluid flow, the thermal resistance of the 
container’s wall, and the effects of natural convection within the 
melt could be ignored for the conditions investigated in this study. 
This model also shows that the rectangular container requires 
nearly half of the melting time as for the cylindrical container with 
the same volume and heat-transfer area. Using the same 
formalism, Najjar and Hasan [73] built a model of a greenhouse 
using a PCM as a heat storage unit. They validate their model by 
comparison with the experimental results of [72]. 

Other studies have also explored the behavior of PCM in 
rectangular geometries [74,75]. They are not explicitly reviewed 
herein but provided for the interested reader. Table 2 presents a 
summary of numerical methods used in the context of the 
Cartesian coordinates system. The symbols FG stand for fixed grid, 
FD for finite difference, FE for finite element, A for analytic, CFD for 
computed fluid dynamics and MM for moving mesh respectively. 


Fujii [72] 


Table 2 
Models for rectangular geometry. 
Model Material Numerical Comment Validation 
formulation 
Shamsundar and FD 2D Surface integrated heat-transfer rates; boundary temperatures 
Sparrow [44,45] solidified fraction and interface position all as function of the time 
Hamdam and n-Octadecane A 1D Propagation and inclination of the interface and energy storage rate Experimental [50] 
Elwerr [49] predicted 
Numerical [214,215] 
Lacroix et al. [51-54] n-Octadecane [51] FG 2D The melting process is chiefly governed by the magnitude of the Experimental [53,246] 
Stefan number 
Numerical [16] 
Brousseau and FD 2D Multi-layer PCM storage Numerical [251] 
Lacroix [55,56] 
FG Parametric model formulated from numerical study 
Costa et al. [57] n-Octadecanol FD 2D For n-octadecanol, there is poor agreement with experimental Experimental 
gallium, tin results. Variation of viscosity with temperature, as heating [245-249] 
losses in the wall or different initial 
supercooling are the proposed explanation for this discrepancy 
Costa et al. [65] Multiple PCM FD 1D For 1D model, calculations have been made for the melt Numerical 
fraction and energy stored for conduction and convection, [216,217] 
while for 2D model, calculations have been made for conduction only. 
2D 
Vakilaltojjar and CaCl2-6H20 A 2D Thin flat containers and air is passed through gaps between them On going, not 
Saman [67] published 
KF-4H20 FE 
Dolado et al. [68] A 1D Encapsulated PCM slab Experimental [218] 
FD 2D 
CFD 
Zukowski [69] Paraffin wax FD 3D Interpolating cubic spline function method is used for determining an Yes 
RII-56 effective specific heat 
Silva et al. [70] Paraffin wax MM 1D Charge and discharge of a PCM encapsulated slab Yes 
Vynnycky and Lauric acid MMLE* 2D Laplace-Euler: new numerical method* Experimental [116] 
Kimura [71] 
Applicable to water and gallium Numerical [219] 
Zivkovic and CaCl2-6H20 FD 1D Isothermal phase change of encapsulated PCM Yes 


Numerical [25] 
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4.2. Spherical geometry 


Spherical geometry represents an interesting case for heat 
storage application, since spheres are often used in packed beds as 
those described in Section 4.4. Due to the complexity of such 
systems, it is often more efficient to first model the behavior of an 
individual sphere, to then describe it with a simple parametric 
model to be used in the packed bed modeling. 

Roy and Sengupta [76] have examined the melting process with 
the solid phase initially uniformly supercooled. The presence of 
supercooling, forced them to modify the heat-transfer equation to 
include the effects of a temperature gradient in the solid core. As a 
result, a closed-form solution could not be obtained. At every time 
step, the unsteady conduction equation has been solved numeri- 
cally using a toroidal coordinate system, which has been suitably 
transformed to immobilize the moving boundary and to transform 
the infinite domain into a finite one. In a subsequent paper [77], 
they have examined the impact of convection on the melting 
process. They demonstrated that the fluid flow in the upper liquid 
region is essentially quasi-steady, since the liquid velocity in both 
the film and the upper zone are much greater than the rate of 
movement of the interface. They also demonstrated the impor- 
tance of the ratio of the buoyancy force due to density variation of 
the liquid with temperature, and the density difference between 
the two phases, which reduces the melt rate and places an upper 
bound for the range over which the film solution is valid. 

Ettouney et al. [78] experimentally evaluated the heat transfer 
in paraffin inside spherical shell. Their calculation shows that the 
Nusselt number in the phase-change material during melting is 
one order of magnitude higher than during solidification and is 
strongly dependent on the sphere diameter. In consequence, the 
melting was about a factor two slower than the freezing. This 
counter-intuitive behavior is caused by the convection, which 
brings a large fraction of the heat to the upper part of the sphere, 
where it heats the melt fraction instead of the solid phase. In the 
freezing process, conduction dominates and the freezing occurs 
around the entire surface. As the melt volume decreases and the 
role of natural convection diminishes rapidly. This result contrasts 
with the work of Barba and Spriga [79] on the behavior of 
encapsulated salt hydrates, used as latent energy storage in a heat- 
transfer system of a domestic hot water tank, who found that the 
shortest time for complete solidification is provided by small 
spherical capsules, with high Jacob numbers and thermal 
conductivities. On the impact of convection, Khodadadi and Zhang 
[80] also noted that it created an asymmetry: melting in the upper 
region of the sphere much faster than in the lower region due to the 
enhancement. Their computational findings were verified through 
qualitative constrained melting experiments using a high-Prandtl 
number wax as the phase-change material. 

Asimple analytical solution was derived by Fomin and Saitoh [81] 
to describe close-contact melting within a spherical capsule. Its 
accuracy is estimated to be 10-15% by comparison to numerical 
analysis. This tool is likely to be useful when evaluating thermal 
energy storage systems which contain thousands of capsules. A 
similar analysis has been done for elliptical capsule. It was found that 
elliptical capsules show higher rate of melting than circular ones. 
Prolate capsules provide more effective melting than oblate ones, 
even though they have the same aspect ratios and vertical cross- 
sectional areas [82,83]. In addition, Adref and Eames [84,85] have 
been able to derive semi-empirical equations that allow the mass of 
ice within a sphere to be predicted during the freezing or melting 
processes. Complementary studies on quasi-static models and on 
parametric models have been done by others research teams [86-90]. 

Ismail and Henriquez [88] created a parametric model for the 
study of the solidification of a PCM in a spherical shell. This model 
was based on pure conduction in the PCM subject to boundary 


conditions of constant temperature or convection heat transfer on 
the external surface of the spherical shell. The model was validated 
by comparison with available similar models [91-94] and the 
agreement was found to be satisfactory. Like other groups, they 
found that the increase of the size of the sphere leads to increasing 
the time for freezing. The same team performed a specific 
numerical and experimental study on water solidification [89]. 
The governing equations of the problem and associated boundary 
conditions were formulated and solved using a finite-difference 
approach and a moving grid scheme. A shortcoming of the model 
was its incapacity of adequately handling supercooling. Neverthe- 
less, they conclude that shell material affects the solidification 
time. Freezing was faster for metallic material (Cu, Al) followed by 
polyethylene, acrylic and PVC. However, the time difference was 
rather small, which still makes non metallic capsules attractive. 

Veerappan et al. [95] did a an investigation the phase-change 
behavior of 65 mol% capric acid and 35 mol% lauric acid, calcium 
chloride hexahydrate, n-octadecane, n-hexadecane, and n-eico- 
sane inside spherical enclosures. They compared their numerical 
model results with the experimental work of Eames and Adref [96] 
and they concluded with a good agreement between them, 
notwithstanding the fact that Eames and Adref [96] worked with 
ice which floats unlike solid phase of modeled PCMs. Based on their 
modeling, they found that spheres with larger diameters take more 
time to completely solidify than smaller diameters ones, since heat 
has to travel a smaller distance and hence rate of solidification is 
higher for smaller capsules. The same is true for melting. However, 
both freezing and melting took approximately the same time, 
which contrasts with result of [69,78]. 

Regin et al. [97] carried out both an experimental and numerical 
analysis of a spherical capsule placed in a convective environment 
filled with wax as a PCM. Model results were compared to 
experimental results and were found to be in good agreement. 
From the model results, they found that the shorter time for 
complete melting and higher time averaged heat flux occurs in the 
capsules of smaller radii and for higher Stefan number. 

Lin and Jiang [86] developed an improved quasi-steady analysis 
model for freezing in a slab, a cylinder, and a sphere. Their 
improvement was the introduction of an additional term to the 
temperature profile to simulate the transient effect on the 
temperature distribution in the solid phase. This additional term 
is based on the ratio of the heat flux at the phase boundary to that 
at the cooling surface, and physically, represents the thermal 
capacity effect in the frozen region. The maximum relative error of 
the moving phase front location obtained from the improved 
quasi-steady analysis of a slab is about 3% in comparison with that 
obtained from the exact solution of the freezing process in a slab, 
while the error for the classical quasi-steady analysis is 20%. Since 
there is no exact solution available for the freezing process taking 
place in a cylinder or a sphere, the results obtained from the 
improved quasi-steady analysis are compared with results of 
[92,98-102] compiled by Lunardini [103]. For these geometries, 
the maximum relative errors of the improved quasi-steady 
analysis are less than 4%, while the maximum relative errors of 
the quasi-steady approximation are higher than 42%. In addition to 
this improvement in precision, the improved model is also less 
sensitive to the variation of Stephan numbers. 

Bilir and Ilken [87] have investigated the solidification problem 
of a phase-change material encapsulated in a cylindrical/spherical 
container. The governing non-dimensional form of equations of the 
problem and boundary conditions were formulated and solved 
numerically by using enthalpy method with control volume 
approach. The numerical code was validated through comparison 
with results of other authors [88,101,104]. In general, there is a 
good agreement. The small differences are due to the omission of 
the thermal resistance of the container and the time dependency of 
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the convection coefficient. Since convection increases with the size of 
the container, differences get larger as the sphere diameter increases. 
Running the simulations for various Stefan and Biot numbers and 
superheat parameter, they produced a regression which enables 
them to determine the coefficient of a function predicting the 
freezing time. The coefficient of correlation of this regression is 0.996, 
which makes it suitable for solving engineering problems. 

Assis et al. [90] produced a parametric investigation of melting 
in spherical shells. Their simulations incorporated such phenome- 
na as convection in the liquid phase, volumetric expansion due to 
melting, sinking of the solid in the liquid, and close-contact 
melting. Their model attempts to solve complete transient 
conservation equations simultaneously for solid PCM, liquid 
PCM, and air, while allowing for PCM expansion, convection in 
the fluid media (melted PCM and air), and solid phase motion in the 
liquid in a similar way as in [105]. This numerical simulation 
compares favorably with results of the experimental investigation 
based on previous work by the same group [106,107]. 

Khodadadi and Zhang [80] performed a computational study of 
the effects of buoyancy-driven convection on constrained melting of 
phase-change materials within spherical containers. They found 
that, early during the melting process, conduction was dominant. As 
the buoyancy-driven convection is strengthened due to the growth 
of the melt zone, melting in the top region of the sphere is much 
faster than in the bottom region due to the enhancement of the 
conduction mode of heat transfer. They noted that the intensity of 
natural convection is more clearly indicated by the Rayleigh number 
than the Stefan number, as the flow pattern and overall convective 
effects change markedly with the Rayleigh number. These numerical 
findings were backed by an experimental visualization of the 
melting of a commercial grade wax inside spherical bulbs. 

Tan et al. [108,109] carried out a similar work on constrained 
melting of phase-change materials to impend the fall of solid phase 
at the bottom of the container. Paraffin wax (n-octadecane) was 
constrained, through the use of thermocouples, during melting 
inside a transparent glass sphere. Their experimental setup provided 
detailed temperature data that were collected during the course of 
the melting process along the vertical diameter of the sphere. They 
pointed out that experimental setup underestimates the waviness 


and the excessive melting of the bottom of the PCM compared to 
numerical predictions. They concluded that this discrepancy could 
be caused by the support structure required to hold the sphere, since 
it could have inhibited heat from reaching the bottom of the sphere. 

Table 3 presents the synthesis for the analyses in spherical 
geometries. 


4.3. Cylindrical geometry 


Anica Trp [110,111] analyzed the transient heat-transfer 
phenomenon during technical grade paraffin melting and solidifi- 
cation in a cylindrical shell. The mathematical model, formulated 
to represent the physical problem, has been proved suitable to 
treat both melting and solidification processes. It can be concluded 
that the selection of the operating conditions and geometric 
parameters dimensions depends on the required heat-transfer rate 
and the time in which the energy has to be stored or delivered. In 
consequence, these parameters must be chosen carefully to 
optimize thermal performances of the storage unit. 

Gong and Mujumdar [112] developed a finite-element model for 
an exchanger consisting of a tube surrounded by an external co-axial 
cylinder made up of PCMs. This model compares characteristics of 
two operation modes. In mode 1, hot and cold fluids are introduced 
from the same end of the tube whereas in mode 2, hot and cold fluids 
are introduced from different ends. Analyses show that the energy 
charged or discharged in a cycle using mode 1 is 5.0% higher than 
when using mode 2. The charge/discharge rate is also faster when 
using mode 1 because the temperature difference between the fluid 
and the PCM is higher in the fluid inlet than in the outlet. The larger 
the temperature difference is the more deeply the phase-change 
interface penetrates into the PCM and the more heat is stocked. In 
the discharge process, symmetrical phenomena occur. 

El-Dessouky and Al-Juwayhel [113] presented a technique to 
predict the effect of different designs and operating parameters on 
the entropy generation number or the second law effectiveness of a 
latent heat thermal energy storage system (LHTES) where the PCM 
is contained in the annulus of a double-pipe heat exchanger. One of 
the conclusions of this analysis is that second law effectiveness 
increases with the increase of Reynolds’s number, specific heat- 


Table 3 
Models for spherical geometry. 
Model Material Numerical Comment Validation 
formulation 
Roy and Sengupta [76] MM 2D Treatment of convection Numerical [220] 
FD Two zones model 
Barba and Spriga [79] 37.5% NH4NO; +62.5% A 1D Transient position of interface, temperature distribution, No 
Mg(NO3)-6H20 melting fraction, energy released, and duration of 
complete solidification 
Fomin and Saitoh [81] n-Octadecane A 2D Contact melting on unfixed solid phase Experimental [221] 
Numerical [222,223] 
Adref and Eames [84,85] Ice Capsule filled a 80% with an air cell on the top Eames and Adref [96] 
Ismail and Henriquez [89] MM FD 1D Numerical correlation relating the working fluid temperature Numerical [91-94] 
to the time has been produced 
Ismail, Henriquez Ice MM FD 1D Analysis of the impact of thermal conductivity of the shell Yes 
and da Silva [88] material 
Veerappan et al. [95] Various PCM A 2D Solidification and melting of sphere with conduction, natural Experimental [96] 
convection, and heat generation 
Regin et al. [97] Paraffin wax FD Convective environment outside capsule Yes 
Lin and Jiang [86] A 1D  Quasi-steady analysis Numerical [92,98-102] 
Bilir and Ilken [87] FG 1D Correlations which express the dimensionless total solidification | Numerical [88,101,104] 
time of the PCM in terms of Stefan Number, Biot Number and 
Superheat Parameter were derived 
Assis et al. [90] Parrafin wax RT27 MM 2D Partly filled capsule with open end Yes + experimental 
[106,107] 
CFD 
Khodadadi and Zhang [80] Beewax FG 2D Constrained melting Yes 
Tan et al. [108] n-Octadecane MM CFD 2D Constrained and unconstrained melting Experimental [109] 
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transfer area and operating fluid wall temperature in the heat 
recovery process and with the decrease of fluid inlet temperature 
in the heat storage process. 

For a solar water heater, Prakash et al. [114] analyzed a storage 
unit containing a layer of PCM at the bottom. Bansal and Buddhi 
[115] did the same thing for a cylindrical storage unit in the closed 
loop with a flat plate collector. In an effort to improve performances 
of phase-change storage units for solar heating applications, Farid 
[116] has suggested using many PCM with different melting 
temperatures. This idea has been applied later by Farid and Kanzawa 
[117] and Farid et al. [118] for a unit consisting of vertical tubes filled 
with three types of waxes having different melting temperatures. 
During heat charge, the air flows first across the PCM with the 
highest melting temperature to ensure continuous melting of most 
of the PCM. The direction of the air flow must be reversed during heat 
discharge. Both the theoretical and experimental measurements 
showed an improvement in performances of the phase-change 
storage units using this type of arrangement. 

Jones et al. [119] performed well-controlled and well-charac- 
terized experimental measurements of melting of a moderate- 
Prandtl-number material (n-eicosane) in a cylindrical enclosure 
heated from the side. The melt front was captured photographi- 
cally and numerically digitized. A numerical comparison exercise 
was undertaken using a multi-block finite-volume method and the 
enthalpy method for a range of Stefan numbers. Very good 
agreement was obtained between the predictions and the 
experiments for Stefan numbers up to 0.1807. 

Esen and Ayhan [120] developed a model to investigate the 
performance ofa cylindrical energy storage tank. In the tank, the PCM 
was packed in cylinders and HTF flowed parallel to it. The PCMs 
considered were calcium chloride hexahydrate (CaCl2—6H20), paraf- 
fin wax (P-116), NazSO4-10H20 and paraffin. The performance of 
cylindrical energy storage tank was determined by computer 
simulations [121,122] and backed by experimental data. The results 
show that the stored energy becomes higher at a given time as the 
mass flowrate or inlet HTF temperature increases and for more energy 
storage appropriate cylinder wall materials and dimensions should be 
selected, such as higher thermal conductivity and small radius. 

Jian-you [123] produced a numerical and experimental 
investigation of a thermal storage unit involving a triplex 
concentric tube with PCM filling in the middle channel, with 
hot heat-transfer fluid flowing into the outer channel during 
charging process and cold heat-transfer fluid flowing into the inner 
channel during discharging process. A simple numerical method, 
according to energy conservation, called temperature and thermal 
resistance iteration method has been developed for the analysis of 
PCM solidification and melting in the triplex concentric tube. 
Comparison between numerical predictions and the experimental 
data shows good agreement. 
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Table 4 presents a complete picture at a glance of numerical 
methods in cylindrical geometries. 


4.4. Packeds beds 


Since most PCMs possess a low thermal conductivity, increasing 
surface/volume ratio is appealing to increase the heat-transfer 
rate. This can be done by packing a volume with a large number of 
PCM capsules. Saitoh and Hirose [124] put in evidence this idea 
both theoretically and experimentally. However, this improve- 
ment comes at the price of a significant pressure drop and a higher 
initial cost. 

Zhang et al. [125] produced a general model for analyzing 
thermal performances of both melting and solidification process of 
a LHTES composed of PCM capsules. This model can be used to 
analyze the instantaneous temperature distribution, the instanta- 
neous heat-transfer rate, and the thermal storage capacity of 
LHTES system. From this modeling effort, it was found that the 
thermal performance of a given system can be determined by the 
effective filling factor of a PCM capsule. Since formulation of the 
effective filling factor for given geometric and flow conditions are 
available in the literature [1,126-128], this largely reduces 
complexity of calculation. This model was validated with the 
results in literature for a LHTES composed of PCM spheres. The 
model is quite general so it can be used to optimize and to simulate 
the thermal behavior of LHTES. 

Benmansour et al. [129] developed a two-dimensional numerical 
model that predicts the transient response of a cylindrical packed 
bed randomly packed with spheres of uniform sizes and encapsu- 
lating paraffin wax with air as a working fluid. Experimental 
measurements of temperature distribution compared favorably with 
the numerical results over a broad range of Reynolds numbers and 
showed that the time at which the bed reaches the thermal 
saturation state decreases significantly with the mass flow rate of air. 

Bédécarrats et al. [130-136] worked on a model of a cylindrical 
tank filled with encapsulated phase-change materials. This model 
had the particularity of examining the impact of supercooling, which 
can seriously deteriorate performances of a LHTES [131,132]. 
Numerical results were validated with an experimental test tank. 
From both numerical and experiment results, it was found that the 
optimum charging mode is obtained in the case of the vertical 
position, where the motions due to the natural convection, are in the 
same direction as the forced convection. Also, they demonstrated 
that an increase in flow rate increases the charge and the discharge 
rates. They also recommend to slightly oversize the tank. 

Cheralathan et al. [137] did also a numerical and experimental 
study of the impact of super cooling and porosity. They came to the 
conclusion that lower inlet HTF temperature reduces the super- 
cooling of the PCM and the total charging time. However, the 


Table 4 

Models for cylindrical geometry. 
Model Material Numerical Comment Validation 

formulation 
Trp [110,111] Parafin wax RT 30 FG 2D Melting + solidification Yes 
Gong and Mujumdar [112] 80.5% LiF + 19.5% CaF FG 2D Co-axial exchanger No 
El-Dessouky and Al-Juwayhel [113] Parafin wax, CaCl2-6H20 A 2D Second law No 
Prakash et al. [114] M 2D Vertical cylinder with PCM at the bottom Yes 
Bansal and Buddhi [115] MM 2D Closed loop with a flat plate collector Yes 
Farid and Kanzawa [117] MM 2D Several MCPs with different melting Yes 
temperatures 
Jones et al. [119] n-Eicosane MM 2D Multi-block FVM Yes 
Esen and Ayhan [120] CaCly-6H20 paraffin, Na2SO4-10H20 FD FG 2D Energy storage reservoir Experimental 
paraffin [224] 

Jian-you [123] n-Hexacosane 1D Triplex concentric tube Yes 


Temperature and thermal resistance 
iteration method 
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reduction of heat-transfer fluid temperature also reduces the 
performance of the refrigeration system. Therefore, there is a trade 
to do between inlet temperature and overall performances. Low 
porosity increase storage capacity of the system at the expense of a 
slower charging/discharging rate. 

Arkar and Medved [138] also studied a cylindrical container 
containing spheres filled with paraffin. Their model was adapted to 
take into account the non-uniformity of porosity and the fluid 
velocity. Both are the consequence of a small tube-to-sphere 
diameter ratio of their system. A comparison of the numerical and 
experimental results confirmed the important role of PCMs 
thermal properties especially during slow running processes. 

Wei et al. [139] explored the impact of geometry (sphere, 
cylinder, plate and tube) of capsules on packed bed performances. 
Effects of the capsule diameter and shell thickness and the void 
fraction were also investigated. They discovered that the heat 
release performance decreases in the order from sphere, to cylinder, 
to plate and finally to tube. Working on another aspect of system 
optimization, Seeniraj and Lakshmi Narasimhan [140] study 
considered the impact of heat leak through sidewalls by convection. 
Their results indicate a significant delay in the solidification process 
with a relatively lower charging rate for the LHTS unit with 
convective heat leak. Based on this result, they derived limiting 
values of heat leak ratio for effective charging of these systems. 

Ismail and Henriquez [141] have created a simplified transient 
one-dimensional model of PCM capsule stocked in a cylindrical tank. 
The convection present in the liquid phase of the PCM is treated by 
using an effective heat conduction coefficient. The solution of the 
differential equations was realized by the finite-difference approxi- 
mation and a moving mesh inside the spherical capsule. The 
geometrical and operational parameters of the system were 
investigated both numerically and experimentally and their influ- 
ence on the charging and discharging times was investigated. Ismail 
and Stuginsky [142] compared four basic groups of models of packed 
beds: the continuous solid phase models, Schumann’s model, the 
single phase models and the thermal diffusion models. Numerical 
tests were carried out to evaluate the time consumed by each of 
them. It was found that these models were essentially equivalent 
except for the computational time. For example, model with a 
thermal gradient inside the solid particles needed more than twenty 
times the computing resources required by Schumann’s model. 

Chen [143,144] also developed a one-dimensional porous- 
medium model. Predictions of this model were compared with 
the lumped capacitance model which assumes a uniform tempera- 
ture in the packed capsules and in the coolant flow and with 
experimental data derived from literature. Results show that the 
one-dimensional model has the advantage of predicting the 
temperature distributions of PCM and coolant. For the study of 
transient response of PCM beds, Goncalves and Probert [145] have 
introduced a new non-dimensional number, the Goncalves number, 
which mean integrated value can be used to characterize the 
thermal response of a thermal energy storing system. 

Adebiyi et al. [146] have developed a model of a packed bed for 
high-temperature thermal energy storage that incorporates several 
features: thermal properties variations with temperature, provi- 
sions for arbitrary initial conditions and time-dependent varying 
fluid inlet temperature, formulation for axial thermal dispersion 
effects in the bed, modeling for intraparticle transient conduction, 
energy storage in the fluid medium, transient conduction in the 
containment vessel wall. Energy extraction can be achieved either 
with flow direction parallel to that in the storage mode (concurrent) 
or in the opposite direction (counter-current). This model also 
allows the computation of the first and second law efficiencies. Yagi 
and Akiyama [147] have also studied the heat-transfer problem at 
high temperatures. Six different materials were tested as PCMs: two 
inorganic and four metallic materials. The metal PCMs were found to 


be excellent for heat storage because of uniform temperature in the 
capsule. Heat-transfer simulation was also conducted for a packed 
bed process of spherical capsules and concurrent flow for heat 
storage and release showed better result for effective use of storage 
heat than counter-current flow. 

Wanatabe et al. [148,149] have explored the impact of 
combining PCMs with various melting points on exergy efficiency 
and charging and discharging rates in a latent heat storage system. 
The heat storage module consisted of horizontal cylindrical 
capsules filled with three types of PCMs. The optimum melting 
point distribution of the PCMs has been estimated from numerical 
simulations and also from simple equations. The fast charging or 
discharging rate leads to high exergy efficiency. Both the 
experimental and numerical results showed some improvement 
in charging and discharging rates by use of multiple PCMs. 

MacPhee and Dincer [150] studied the melting and freezing of 
water in spherical capsules using ANSYS, GAMBIT, and FLUENT 6.0. 
These models were validated through a rigorous time and grid 
independence tests and their numerical results were successfully 
compared to the experimental ones of [78]. Energy efficiency was 
calculated and found to be over 99%, though viscous dissipation was 
included. Using exergy analysis, the exergetic efficiencies are 
determined to be about 75% to over 92%, depending on the HTF 
scenario. The two efficiencies decreased when hot transfer fluid flow 
was increased, due mainly to increasing heat losses and exergy 
dissipation. These results are different from those obtained by 
Wanatabe et al. [148,149]. Higher temperature differences between 
hot fluid and melting temperature of PCMs were found to be most 
optimal exergetically, but least optimal energetically due to entropy 
generation accompanying heat transfer. Higher temperature 
differences also increased the rates of charging/discharging. 
However, increasing the HTF flow rate did not decrease the charging 
and discharging times by any relatively significant amount. This 
result is somewhat discrepant with various previous studies [129- 
136]. The same team also developed a simple analytical model [151] 
that describes heat transfer, fluid flow, and thermodynamic, which 
brought essentially the same conclusions as the numerical model. 

A summary for packed beds models is provided in Table 5. 


4.5. Finned geometry 


Another approach to increase the heat exchange efficiency in 
PCMs is to use fins in the heat exchange system. These somehow 
more complex geometries have been extensively studied. 

Sasaguchi et al. [152,153] studied experimentally and theoreti- 
cally effects of the configuration of a finned tube on the heat-transfer 
characteristics of a latent heat thermal energy storage system. They 
concluded that performance of the unit is almost the same for any 
configuration with the same surface area even if it is quite different. 

Lacroix [154] developed a theoretical model, using an enthalpy- 
based method coupled to a convective heat transfer of a shell-tube 
thermal storage unit with PCM on the shell side and the HTF 
circulating inside the tubes. The numerical model was validated 
through comparison with experimental data. In addition, compar- 
ison was done between models of bare and finned tubes. Results 
showed that the annular fins are the most effective for moderate 
mass flow rates and small inlet temperatures. Lacroix and 
Benmadda [155] also studied the behavior of a vertical rectangular 
cavity filled with PCM. They used a fixed-grid model which was 
validated with experimental data. They found that both solidifica- 
tion and melting rates were improved by long fins. However, 
melting rates reach a maximum when the number of fins is 
increased since they interfere with the convection. 

Zhang and Faghri [156] studied the heat-transfer enhancement 
produced by externally finned tube in a latent heat energy storage 
system. The heat conduction in the tube wall, fins and the PCM 
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Table 5 
Models for packed bed. 
Model Material Numerical Comment Validation 
formulation 
Zhang et al. [125] A 1D A general model, great reference paper Numerical [225-229] 
Benmansour et al. [129] Parafin wax FD 2D Mass flow rate is determinant Yes 
Bédécarrats et al. [130-136] PM 2D Supercooling studied Experimental [131,135,230] 
Cheralathan et al. [137] Ice 1D  Supercooling is reduced by lower inlet HTF temperature Yes 
Arkar and Medved [138] Parafin wax RT20 PMFDFG 2D Continuous solid phase model Yes 
Wei et al. [139] Paraffin wax FNP-0090 FD FG 1D Geometry analysis of capsules Yes 
Seeniraj and Lakshmi Heat leaks through sidewalls Experimental [144] 
Narasimhan [140] 
Ismail and Henriquez [141] Ice FD MM 1D Effective heat conduction coefficient Yes 
Ismail and Stuginsky [142] 1D Parametric model 
Chen [143,144] Ice 1D Solution by Laplace transform Experimental [144] 
Goncalves and Probert [145]  MgCl2-6H20 Their own dimensionless number! Yes 
Adebiyi et al. [146] ZrOz First and second law efficiency, high temp Yes 
Cu 
Yagi and Akiyama [147] KNO3-NaNO3 1D Six materials with various melting points Yes 
NaCl, Pb 
Al-Si 
Wanatabe et al. [148,149] 1D Second Law model Yes 
MacPhee and Dincer [150] 2D First and second law efficiency. Fluent Experimental [78] 


were described by a temperature transforming model using a 
fixed-grid numerical method [27]. This model assumes that the 
melting process occurs over a range of phase-change tempera- 
tures, but it also can be successfully used to simulate the melting 
process occurring at a single temperature by taking a small range of 
phase-change temperature. This model was validated through 
comparison with Lacroix [154] and agreement between the two 
was found to be quite good. The results show that the local Nusselt 
number inside the tube cannot be simply given by the Graetz 
solution. The model also shows that transverse fins are a more 
efficient way to enhance the melting heat transfer if the initial 
subcooling exists in the PCM. The height of the fins has significant 
effect on the melting front on the two sides of the fins but has no 
significant effect on the melting front between the transverse fins. 

Velraj et al. [157] studied the impact of internal longitudinal 
fins on a cylindrical vertical tube filled with paraffin. A theoretical 
model that accounts for the circumferential heat flow through the 
tube wall was developed using enthalpy formulation and was 
employed in conjunction with the fully implicit finite-difference 
method to solve the solidification in the convectively cooled 
vertical tube. They also generalized the enthalpy formulation of 
Date [158] to the case where the material has a range of phase- 
change temperature. This model was validated with experimental 
data. From this analysis, they concluded that the solidification time 
is approximately 1/n time the no fin case, where n is the number of 
fins. They also pointed out that, in the central part, the flow is 
restricted between fins. Therefore, when the number of fins 
exceeds 4, half the fins can extend to only half the cylinder radius. 
In their subsequent work [159], they pointed out the necessity to 
include the effect of circumferential heat flow through the tube 
wall for higher values of Biot numbers in order to correctly predict 
the heat-transfer behavior. For lower Biot numbers, addition of fins 
makes the surface heat flux more uniform, whereas for higher Biot 
numbers the addition of fins improves the magnitude of the 
surface heat flux and appreciably reduces the solidification time. 
For a given quantity of heat to be extracted, a combination of lower 
Biot number and higher Stefan number is recommended for the 
uniform extraction of heat. The same team [160] also investigated 
transient behavior of a finned-tube latent heat thermal storage 
module with thin circumferential fins. Numerical results indicate 
an appreciable enhancement in the energy storage process with 
the addition of fins in the module. 

Lamberg and Sirén [161-164] developed an analytical model 
for analysis of the melting and solidification in a finned PCM. The 


results of the derived analytical model were compared to 
numerical results of a two-dimensional model solved numerically 
by means of the effective heat capacity method and enthalpy 
method. The two-dimensional results were compared to simplified 
one-dimensional results in order to find out the accuracy of the 
simplified analytical model. The results of the experimental work 
were then compared to the numerical results calculated using the 
effective heat capacity method and enthalpy method, in order to 
find out the accuracy of different numerical methods. The results 
show that the analytical models give a satisfactory estimate of fin 
temperature and of the solid-liquid interface. 

Castell et al. [165] studied both experimentally and numerically 
the impact of external longitudinal fin on PCM module placed in a 
water tank working witha solar heating system. The geometry of the 
PCM modules is the same used in an other experimental work done 
by the authors [166,167] but with added fin on the PCM module. 
Experimental work was carried out to determine natural convection 
heat-transfer coefficients for PCM modules. To compare the 
behavior of the system using finned PCM modules with another 
one without fins, two numerical codes were implemented. The 
results obtained from the simulation were compared with the 
experimental ones to validate the numerical codes. Based on those 
experimental results, useful Nusselt correlations in function of 
Rayleigh number were found to evaluate the natural convection 
heat-transfer coefficient for that specific geometry. 

Reddy [168] modeled a PCM-water solar system consisting of a 
double rectangular cross section enclosure where the top enclosure 
is filled with paraffin wax and the bottom is filled with water. The 
numerical modeling has been carried out using the commercial CFD 
software FLUENT. The geometry modeling and mesh were generated 
in GAMBIT. An enthalpy porosity technique was used in FLUENT for 
modeling the solidification/melting process. In this technique, the 
mushy zone is modeled as a pseudo porous medium in which the 
porosity decreases from 1 to 0 as the material solidifies. Simulation 
results indicate that an optimal number of fins should be used to 
optimize the overall system performance. 

Gharebaghi and Sezai [169] investigated the enhancement of 
energy storage rate of a thermal energy storage unit filled with a 
PCM by inserting a fin array. Heat is transferred to the unit through 
the container walls, to which fins are attached; a commercial 
paraffin wax is stored between the fins. A mathematical model, 
based on the finite-volume method for a two-dimensional domain, 
is developed for solving the melting problem. A fixed-grid method 
is used to simulate the melting of the PCM. The interface between 
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the solid and the melt is modeled as a porous medium. Transient 
simulations have been performed on non-uniform grids for 
different fin and PCM layer thicknesses. For horizontal and vertical 
arrangements of modules, it was observed that for high tempera- 
ture differences, the heat-transfer rate can be increased as much as 
80 times by adding a fin array. The minimum enhancement rate 
was 3 times, which was achieved with widely spaced fins. It was 
also observed that the Nusselt number is higher for vertical 
modules arrangement compared with that of horizontal modules 
for all fin spacings and temperature differences. 

Ismail et al. [170] studied a thermal numerical model for the 
solidification of PCM around a radially finned tube with a constant 
wall temperature. Their model was based on a pure conduction 
formulation and the enthalpy method. Numerical simulations 
were performed to investigate the effects of the number of fins, fin 
thickness, fin material, aspect ratio of the tube arrangement and 
the tube wall temperature. The study shows that the geometrical 
aspects have strong influences on the time for complete solidifica- 
tion and on the solidification rate. Also, the tube material as well as 
the tube wall temperature affect the time for complete solidifica- 
tion. The same team [171], following the same procedure, analyzed 
the behavior of the solidification of a PCM around a vertical axially 
finned isothermal cylinder. Their numerical model was validated 
by experimental data. From this analysis, the authors suggested 
that a metallic tube fitted with four or five fins of constant 
thickness equal to the tube wall thickness and of radial length 
around twice the tube diameter should be the best compromise 
solution between efficiency, increase in the heat flow rate and the 
loss of available storage capacity. 

Kayansayan and Ali Acar [172] modeled the formation of ice 
around a finned-tube heat exchanger. Due to the small tempera- 
ture difference between the tube outside surface and the PCM (less 
than 5°C), natural convection is very weak and even more 
hampered by the presence of vertical fins. The validity of their 
model was established by comparison with Lacroix [154] and an 
experimental setup. 

The thermal control of portable electronic devices, with their 
intermittent operation, is well adapted to PCM utilization. 
However, careful optimization needs to be done to insure the 
adequate sizing and consistent behavior under realistic operation- 
al condition. 

Akhilesh et al. [173] developed a model of a composite heat sink 
composed of vertical arrays of fins surrounded by phase-change 
material. It should be noted that heat load was coming from the top 
of the heat sink. This heat sink was to be used for cooling electronic 
package below a set point temperature. This model provides a 
thermal design procedure for proper sizing of the heat sink, for 
maximizing energy storage and for operation time. Based on a 
scaling analysis of the governing two-dimensional unsteady 
energy equations, a relation between the critical dimension for 
the heat sink and the amount of PCM used was developed. 

Shatikian et al. [105,174,175] explored numerically the process 
of melting of a phase-change material (PCM) in a heat storage unit 
with internal fins open to air at its top. The phase-change material, 
paraffin wax, was stored between the fins. Transient three and 
two-dimensional simulations were performed using FLUENT 6.0, 
yielding temperature evolution in the fins and the PCM. Using a 
dimensional analysis to generalize their results, they have shown 
that within the same geometry, the Nusselt number and melt 
fraction depend on the product of the Fourier and Stefan numbers. 
For relatively wide vertical PCM layers, the Rayleigh number must 
also be included, to take into account the effects of convection at 
advanced stages of the melting process. 

Béhunek and Fiala [176] also modeled a heat sink for electronic 
package including PCM based on the work of Nayak et al. [177]. In 
their configuration, the inorganic crystalline salt acting as PCM was 


placed in a chamber right over the chip. This chamber was put in 
contact with fin cooled by forced convection of air. They first 
created an analytical description and a solution of heat transfer, 
melting and freezing process in 1D and compared their solution to 
the numerical solution obtained, through a FEM solution in ANSYS 
software, of a real 3D cooler. In addition, results of 3D numerical 
solutions were verified experimentally. 

Saha et al. [178] also made an analysis of a similar system and 
geometry both experimentally and numerically. Their experimen- 
tal setup consisted of an electrical heater to simulate the electronic 
chips with the heat sink filled with eicosane above. The governing 
equations of their numerical model follow a single domain 
approach where both fluid flow and heat-transfer equations are 
solved simultaneously. The fin, heater, and substrate materials are 
assigned very high viscosity and melting points, which ensures 
that they remain solid. In consequence, a common set of governing 
conservation equations for both the solid and liquid regions can be 
used [177]. An optimal fin fraction of 8% was found. This was 
attributed to convection cells in the molten region of the PCM 
present at a large volume fraction of PCM but otherwise 
suppressed. It was also found that a large number of narrow 
cross section fins is preferable, since it improves the thermal 
contact for a same volume. 

Kandasamy et al. [179] have also studied a similar heat sink 
experimentally and numerically. Results show that increased 
power inputs enhance the melting rate as well as the thermal 
performance of the PCM-based heat sinks until the PCM is fully 
melted. A three-dimensional computational fluid dynamics model 
was used to simulate the problem. To describe the PCM-air system 
with a moving internal interface but without inter-penetration of 
the two fluids, the volume-of-fluid model has been used 
successfully [180] and showed good agreement with experimental 
data. The same team [181] using the same approach explored the 
effect of orientation on performances. They concluded that 
orientation has a minimal effect in their configuration. 

In an original approach, Wang et al. [182] created a constructal 
rule for the design of fins in PCM application. Constructal 
optimization of heat conduction with phase-change mimics the 
growth of plants in nature. The high conductivity material is 
treated as the plant root, and the PCM acts as the soil. The product 
of thermal conductivity and temperature gradient integration 
over time in the whole melting time interval of the PCM is taken as 
the criterion to determine the arrangement position of the high 
conductivity material. Its physical meaning is the quantity of total 
heat transported during the melting time interval. The generation 
rule of constructal optimization of heat conduction with phase 
change is to grow the high conductivity material at the position 
with the maximum total heat transported during the time 
interval. The degeneration rule is to withdraw the high 
conductivity material at the position with minimum total heat 
transported during the time interval. Comparison between 
constructural and man-made shapes confirm that constructal 
ones are superior in melting time and in the average cold release 
power. 

Table 6 presents the synthesis for finned surfaces. 


4.6. Porous and fibrous materials 


Erk and Dudukovic [183] modeled and experimented with a 
PCM consisting of n-octadecane retained by capillary forces in a 
porous silica support. It optically changes from opaque to 
semitransparent when the wax melts, thereby allowing the melt 
front within the bed to be tracked. This property and the outlet 
temperature were compared with the model. Measured outlet 
temperature compares qualitatively with model predictions. The 
model quantitatively predicts melt front movement in the first 60% 
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Table 6 
Models for finned surfaces. 
Model Material Numerical Comment Validation 
formulation 
Sasaguchi et al. [152,153] 3D Performance independent of geometry only surface area Yes 
Lacroix [154] n-Octadecane FD 2D An analytic 1D model is also produce Yes 
Lacroix and Benmadda [155] FG 
Zhang and Faghri [156] FD 2D Validation of semi-analytic result of [231] Numerical Lacroix [154] 
FG Kays and Crawford [232] 
Velraj et al. [157,159,160] Paraffin RT 60 [157] FD 2D Phase-change material kept inside a longitudinal Yes [157] 
LiF-MgF, [160] FG internally finned vertical tube 159] num. [155] 
[160] num. [233] 
Lamberg and Sirén [161-164] Paraffin FE 2D Analytical model validated Yes 
Castell et al. [165,234] Sodium acetate A Horizontal and vertical fins around circular vertical cylinder Yes 
trihydrate 
with graphite 
Reddy [168] Paraffin wax MM 2D Double rectangular enclosure where the top enclosure is Experimental [252] 
filled with paraffin wax and the bottom is filled with water 
Gharebaghi and Sezai [169] Paraffin RT27 FD 2D Slabs containing paraffin and metal fins made of aluminum Yes 
FG 
Ismail et al. [170,171] Ice [105] FD 2D Radially finned horizontal tube [170] Yes [171] 
Paraffin [174] Axially finned vertical tube [171] 
Kayansayan and Acar [172] Ice FD 2D Radially finned horizontal tube Numerical Lacroix [154] 
FG Zhang and Faghri [156] 
Akhilesh et al. [173] n-Eicosane FD 2D Array of vertical fins Numerical Bejan [235] 
FG 
Shatikian et al. [105,174,175] Paraffin wax RT25 2D Array of vertical fins Numerical [105] 
3D 
Béhunek and Fiala [176,163] CaClz-6H20 1D 3D Also analytical 1D Yes 


of the bed. Discrepancies between the model and experiments are 
linked to significant heat losses in the small insulated bench-scale 
apparatus used for the experiment.Mesalhy et al. [184,185] studied 
numerically and experimentally the effect of porosity and thermal 
properties of a porous medium filled with PCM. In their model, the 
governing partial differential equations describing the melting of 
phase-change material inside porous matrix were obtained from 
volume averaging of the main conservation equations of mass, 
momentum, and energy. Due to the huge difference in thermal 
properties between the phase-change material and the solid 
matrix, two energy equations model was adopted to solve the 
energy conservation laws. This model can handle local thermal 
non-equilibrium condition between the PCM and the solid matrix. 
Finite-volume technique in conjunction with numerical grid 
generation based on body fitted coordinate transformation has 
been adopted. From their model, it was found that the best 
technique to enhance the response of the PCM storage is to use a 
solid matrix with high porosity and high thermal conductivity. The 
same team [186] investigated the properties of graphite foams 
infiltrated with phase-change materials. The geometry studied was 
two concentric cylinders, the inner being kept at a constant 
temperature while the outer one is filled with foam/PCM matrix. 
The model has been validated by applying it on the case described 
by Khillarkar et al. [187] and has shown to be in good agreement. 
Model results indicate that estimated value of the average output 
power using carbon foam of porosity 97% is about five times 
greater than that for using pure PCMs. 

Fukai et al. [188,189] modeled and experimented the use of 
carbon-fiber combined with n-octadecane. Carbon fibers were 
added to improve the thermal conductivities of phase-change 
materials packed around heat-transfer tubes. The transient 
thermal responses improved as the volume fraction of the fibers 
and the brush diameter increased. However, these do not further 
improve when the diameter of the brush is larger than the distance 
between the tubes due to thermal resistance near the tube walls. 
The discharge rate using the brushes with a volume of one percent 
was about 30% higher than that using no fibers. In the charge 
process, the brushes prevent natural convection. Still, the charge 


rate with the brushes is 10-20% higher than that with no fibers. A 
heat-transfer model describing anisotropic heat flow in the 
composite was numerically solved. The calculated transient 
temperatures agreed well with the experimental data. A simple 
model was also developed to predict the heat exchange rate 
between the composite and the heat-transfer fluid. Hamada et al. 
[190] compared experimentally and numerically the thermal 
behavior of carbon chip and carbon brush. The carbon-fiber chips 
are efficient for improving the heat-transfer rate in PCMs. 
However, the thermal resistance near the heat-transfer surface 
is higher than that for the carbon brushes. As a result, the overall 
heat-transfer rate for the fiber chips is lower than that for the 
carbon brushes. Consequently, the carbon brushes are superior to 
the fiber chips for the thermal conductivity enhancement under 
these experimental conditions. 

Elgafy and Lafdi [191] studied experimentally and analytically 
the behavior of nanocomposite carbon nanofibers filled with 
paraffin wax. An analytical model was introduced based on a one- 
dimensional heat conduction approach to predict the effective 
thermal conductivity for the new nanocomposites. Their findings 
showed good agreement with the experimental data of the authors 
and of Lafdi and Matzek [192]. 

Frusteri et al. [193] studied experimentally and numerically 
thermal conductivity and charging-discharging kinetics of 
inorganic PCM44 containing carbon fibers. A numerical model, 
based on the enthalpy formulation and on an implicit finite- 
difference method, has been developed to predict the one- 
dimensional heat-transfer diffusion of a PCM containing carbon 
fibers. The results obtained by use of this model were satisfactory 
compared to experimental results. Authors also pointed out that 
the values of thermal conductivity calculated by Fukai’s method 
are significantly different from those measured experimentally. 
This was probably caused by a reduction of the thermal 
conductivity of the composite material near the heat-transfer 
surface due to both the wall thermal resistance and the contact 
thermal resistance [190]. 

The synthesis for porous and fibrous materials is given in 
Table 7. 
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Table 7 
Models for porous and fibrous materials. 
Model Material Numerical Comment Validation 
formulation 
Erk and Dudukovic [183] n-Octadecane MM 1D Second law efficiency calculated Yes 
Mesalhy et al. [184,185] FD 2D Carbon foam impregnated with PCM Khillarkar et al. [187] 
FG 
Fukai et al. [188,189] n-octadecane FD 2D Carbon-fiber brushes in PCM Yes 
FG 
Hamada et al. [190] n-octadecane FD 2D Carbon-fiber chips in PCM Yes 
FG 
Elgafy and Lafdi [191] Parafin wax A 1D Carbon nanofibers in paraffin wax Lafdi and Matzek [192,191] 
Frusteri et al. [193] 37:26:38 2D Carbon fibers in PCM Yes 
NH,NO3 
Mg(NO3)2‘6H20 
MgCl2-6H20 


4.7. Slurry 


Like porous and fibrous material, slurries increase heat-transfer 
performance and capacity of fluid. This enhancement improves the 
overall efficiency of cooling or heating system. Reader should 
consult the review of Zhang et al. [194] to learn more about 
properties and applications of phase-change material slurries. 

This research field was first explored by Charunyakor et al. 
[195]. They formulated a model for heat transfer of microencap- 
sulated phase-change material slurry flow in circular ducts. Heat 
generation or absorption due to phase change in the particles was 
included in the energy equation as a source term. The enhance- 
ment of thermal conductivity due to the particle/fluid interactions 
was also taken into consideration. Results of this model have 
demonstrated that Stefan number and PCM microcapsule concen- 
tration are the most important parameters for the system 
performances. 

Zhang and Faghri [196] studied laminar forced convection heat 
transfer of a microencapsulated phase-change material in a 
circular tube with constant heat flux. Melting in the microcapsule 
was solved by a temperature transforming model instead of a 
quasi-steady model. This model was validated with results of 
Charunyakorn et al. [195] and Goel et al. [197]. This analysis 
pointed out that the most significant reason for the large difference 
between Goel et al. [197] experimental results and Charunyakorn 
et al. [195] numerical results is that melting takes place over a 
range of temperatures below the melting point. It also put into 
evidence the importance of shell heat-transfer properties to 
explain this discrepancy. 

Alisetti et al. [198,199] introduced an effective heat capacity 
model for heat transfer in PCM slurries. Their objective was to 
address the shortcoming of previous models using complicated 
source terms or special analytical techniques or variation of the 
specific heat of material with temperature. This new formulation 
was easier to implement in standard computer fluid dynamic 
packages. 

Roy and Avanic [200] developed a turbulent heat-transfer 
model based on the effective specific heat capacity for PCM slurry 
in a circular tube with constant wall heat flux. Their model was 
compared with previously published numerical [237,238] and 
experimental data [239-241]. It demonstrated that the bulk Stefan 
number, the melt temperature range, and the degree of subcooling 
are the driving parameters of the system. 

Hu and Zhang [201] presented an analysis of the forced 
convective heat-transfer enhancement of microencapsulated 
phase-change material slurries flowing through a circular tube 
with constant heat flux. The model used an effective specific heat 
capacity and was validated with the results of Goel et al. [197]. It 
was found that the conventional Nusselt number correlations for 
internal flow of single phase fluids are not suitable to accurately 


describe the heat-transfer enhancement with microencapsulated 
phase-change material suspensions, and a modification is intro- 
duced that enables evaluation for the convective heat transfer of 
internal flows. Results from the model show that Stephan number 
and volumetric concentration of microcapsules are the most 
important parameters governing the heat-transfer enhancement 
of phase-change slurries. However, initial supercooling, phase- 
change temperature range and microcapsule diameter also 
influence the degree of heat-transfer enhancement. The enhance- 
ment increases with decreasing initial supercooling and/or phase- 
change temperature range, and with microcapsule diameter. 
Therefore, the degree of enhancement in a small-diameter tube 
may be better than in a large pipe for a given phase-change slurry. 

Hao and Tao [202] developed a model for simulation of the 
laminar hydrodynamic and heat-transfer characteristics of sus- 
pension flow with micro-nano-size PCM particles in a micro- 
channel. The model solves simultaneously the two-phase equa- 
tions of mass, momentum, and energy. This has allowed the 
investigation of heat-transfer enhancement limits, due to non- 
thermal equilibrium melting of PCM, existence of a particle- 
depleted layer, particle-wall interaction, phase-change region 
propagation, and local heat-transfer coefficients for both constant 
wall heat flux and constant wall temperature boundary conditions. 
The same group [203] developed a two-phase, non thermal 
equilibrium-based model for the numerical simulation of laminar 
flow and heat-transfer characteristics of PCM slurry in a micro- 
channel. The model was validated against the experimental results 
of Goel et al. [197]. The results show that for a given Reynolds 
number, there exists an optimal heat flux under which the 
effectiveness is the highest. 

Ho et al. [204] examined the effects of wall conduction on the 
heat-transfer characteristics of solid—liquid phase-change material 
suspension flow. A mixture continuum approach was adopted in 
the formulation of the energy equation, with an enthalpy model 
describing the phase-change process in the phase-change material 
particles. From numerical simulations via the finite-volume 
approach, it was found that the conduction heat-transfer 
propagating along the pipe wall resulted in significant preheating 
of the suspension flow in the region upstream of the heated 
section, where melting of the particles may occur and therefore the 
contribution of the latent heat transfer to convection heat 
dissipation over the heated section was markedly attenuated. 

Inaba et al. [205] studied the fluid flow and heat-transfer 
characteristics for Rayleigh-Bénard natural convection of non- 
Newtonian phase-change material. From their model, they found 
that Rayleigh number, Prandtl number and aspect ratio were the 
most important factors for evaluating the natural convection in 
enclosures for most of Newtonian and non-Newtonian fluids. 
However, some modifications are necessary for evaluating the 
natural convection in a PCM slurry. In consequence, a modified 
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Table 8 
Models for Slurries. 
Model Material Numerical Comment Validation 
formulation 
Charunyakor et al. [195] FD 1D Circular ducts Ahuja [236] 
Zhang and Faghri [196] FD FG 1D Circular tubes, constant heat flux Experimental 
Goel et al. [197] 
Alisetti et al. [198,199] CFD Effective heat capacity 
Roy and Avanic [200] 10% n-hexadecane FD 3D Turbulent heat transfer Numerical Petukhov [237] 
in water FG Sparrow [238] 
Experimental 
Hartnett [239] 
Choi [240,241] 
Hu et al. [201,226] FD 2D Circular tubes, constant heat flux Experimental 
FG Goel et al. [197] 
Ho et al. [204] FD Micro and nano particles Experimental [197] 
FG Numerical [195,196] 
Inaba et al. [205,206] 30:5:65 FD 2D Natural convection, Newtonian and Yes 
Parafin surfactant FG non-Newtonian fluids Numerical 
water Vahl Davis and Jones [242] 
Ozoe and Churchill [243] 
Cassidy [207] Octacosane FD 2D Mixed convection Yes 
FG 
Royon and Guiffant [208] FD 1D Millimetric particles in oil Numerical Bird et al. [244] 
FG 
Sabbah et al. [209] CFD 3D 3D single phase flow, FLUENT Experimental 
Goel et al. [197] 
Kuravi et al. [210] FE 3D Bulk fluid with adaptative specific heat Experimental 
Goel et al. [197] 
Zhang et al. [211,212] FD 2D Analogy with thermal sources Charunyakorn et al. [195] 
FG Alisetti and Roy [199] 
Sun et al. [213] Paraffin wax 2D Neutrally buoyant ceramics Experimental 


Jones et al. [118] 


Stefan number was defined to address this issue. The same team 
[206] produced a two-dimensional numerical simulation of 
natural convection in a rectangular enclosure with non-Newtonian 
PCM microcapsulate slurry. The formulation of the numerical 
model has been done using the finite-volume method. This study 
has demonstrated that heat-transfer coefficients of the PCM slurry 
with phase change are higher than those without phase change. It 
was found that the natural convection effect is intensified, and the 
heat transfer is enhanced in the case of the PCM slurry, due to the 
contribution of the latent heat transfer. 

Cassidy [207] studied numerically and experimentally a micro 
PCM fluid in a circular tube with mixed convection. The PCM was 
octacosane encapsulated by a polyethylene shell to form a 
spherical particle with an average diameter of 20 wm. The micro 
PCM particles were suspended in a 50/50 ethylene glycol water 
mixture. Experimental measurements were made of the tube outer 
wall at the top and bottom of the copper tube. Numerically an 
incompressible flow model was used with an Eulerian-Eulerian 
method to solve the two-phase momentum and energy equations. 
The numerical model was validated using experimental data of 
single phase flow with mixed convection available in the literature 
and was also verified and thermal results of both single phase and 
two-phase flow from the experimental work of the same author. 

Royon and Guiffant [208] constructed a model describing the 
thermal behavior of slurry of phase-change material flow with 
millimetric particles dispersed in oil in a circular duct. They 
formulated a phenomenological equation which takes into account 
the heat generation due to phase change in the particles. A plateau 
of the temperature appeared in all of the simulations where the 
PCM particles were considered. The length of this plateau is 
strongly dependent on both the weight fraction of the particles and 
the duct wall temperature. 

Sabbah et al. [209] has investigated the influence of micro- 
encapsulated phase-change material on the thermal and hydraulic 
performance of micro-channel heat sinks. They developed a three- 


dimensional, one-phase, laminar flow model of a rectangular 
channel using water slurry of MEPCM with temperature depen- 
dent physical properties. The conservation equation was solved 
using FLUENT 6.2. The numerical results have been compared to 
experimental data of Goel et al. [197] once adapted for the circular 
geometry of this experiment. The model results showed a 
significant increase in the heat-transfer coefficient that is mainly 
dependant on the channel inlet and outlet temperatures and the 
selected PCM melting temperature. The enhancement is higher for 
low concentration of PCM in slurry due to the increase of slurry 
viscosity with its concentration. 

Kuravi et al. [210] produced a numerical investigation of the 
performance of a nano-encapsulated phase-change material slurry 
in a manifold micro-channel heat sink. The slurry was modeled as a 
bulk fluid with varying specific heat. The temperature field inside 
the channel wall is solved three dimensionally and is coupled with 
the three-dimensional velocity and temperature fields of the fluid. 

Zhang et al. [211] analyzed the convective heat-transfer 
enhancement mechanism of microencapsulated PCM slurries 
based on the analogy between convective heat transfer and 
thermal conduction with thermal sources. The heat-transfer 
enhancement for laminar flow in a circular tube with constant 
wall temperature is analyzed using an effective specific heat 
capacity model. This model was validated by comparison with the 
result of Charunyakorn et al. [195] and Alisetti and Roy [199]. The 
numerical simulation results have pointed out the Stephan 
number and the phase-change temperature range as the most 
important parameters influencing the Nusselt number fluctuation 
profile of phase-change slurry. The fluctuation amplitude increases 
with decreasing phase-change temperature range. The fluctuation 
amplitude and range increase with decreasing Stephan number 
when phase-change happens. Later, members of the same team 
[212] studied experimentally and numerically characteristics of 
microencapsulated PCM flowing in a circular tube. The enhanced 
convective heat-transfer mechanism was analyzed by using the 
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enthalpy model. Three kinds of fluid-pure water, micro-particle 
slurry and MPCMs were numerically investigated. They had 
essentially reached the same conclusions. 

Using the same setup as Jones et al. [119], Sun et al. [213] 
studied the behavior of slurry composed of neutrally buoyant 
ceramic hollow spheres suspended in a paraffin wax. The 
numerical analysis employed a particle-diffusive model and the 
enthalpy method. Reasonable agreement was obtained between 
the experiments and numerical predictions. Nevertheless, flow and 
heat-transfer characteristics were greatly altered due to the 
presence of the solid particles and that the particle-diffusive 
model was insufficient to describe the particle migration during 
melting. 

A synthesis of the numerical study of slurries is provided in 
Table 8. 


5. Conclusion 
5.1. Content 


This paper presented an exhaustive review of numerical 
methods applied to the solutions of heat-transfer problems 
involving phase-change materials for thermal energy storage. 
The review is a model collection of fundamental and most recent 
works published on the subject. This survey is organized according 
to the problem geometry (Cartesian, spherical, and cylindrical) and 
specific configurations or applications (packed beds, finned 
surfaces, porous and fibrous materials, slurries). The authors do 
not claim anything about the completeness of the review as some 
papers may have been unfortunately neglected. The authors 
apology for any crucial omission and welcome all authors to 
contact them to complete the review. 


5.2. Analysis: numerical issues and validation 


As a general conclusion, one can readily observe that the older 
studies all had experimental counterparts to validate the modeling 
of the problem with an appropriate set of equations and the 
subsequent formulation of a numerical method to solve the 
relevant sets of discretized equations. Many of these early studies 
also involved analytical solutions used to validate the model for 
selected problems that admit closed-form solutions. 

As time passes, the authors rely more and more on other 
studies, mostly other numerical studies, to validate their own 
numerical results. And it is somewhat interesting to mention that 
among the 250+ references cited here, in only one the authors 
stated that the results were not “in good agreement with those 
found in the references”. Many of the recent studies discuss their 
results qualitatively only as the comparison with a graph taken 
from a publication may be somehow hazardous. 

In recent studies, the proportion of analyses which rely on 
commercial codes increases and the discussions that pertain to 
stability, convergence, grid independence and other related 
numerical issues diminishes. 

On the other hand, we found that some laboratories are involved 
in LHTES for decades and in these labs there is a mature culture of 
experimental, analytical, and numerical methods for the analysis. 


5.3. Analysis: type of materials and configurations 


We would like to warn readers about the risk involved when 
extrapolating from numerical and experimental studies and 
making general conclusion from them. For example, convection 
is strongly sensitive to the geometry and the size of the enclosure 
as to the viscosity and thermal conductivity of PCM. Thermal 
conductivity of the enclosure and the PCM solid-liquid density and 


conductivity contrasts also influence the thermal behavior in 
important and counter-intuitive ways. 

There is an incentive to mix several materials with different 
phase-change temperature to allow for more heat recovery. Packed 
beds, finned geometries, porous and fibrous material-based 
applications, and slurries are the particular configurations for 
most studies that pertain to systems involving PCMs. The latter 
receives more attention these days as a result of the incentive to 
study nano-heat transfer in OECD countries. 

However, it would have been interesting in these studies to 
discuss the issue of pumping power or the effect of the micro- or 
nano-capsules of PCMs on the apparent or effective viscosity of the 
fluids. None of the authors mention this critical issue when it 
comes to establish the net energy balance of operation. 

Moreover, in several studies, carbon fibers are used to enhance 
conduction in PCMs (it is a well known fact that their average 
conductivity is poor) but these studies are silent about the energy 
cost of production of those fibers. 

Another important aspect often neglected when modeling 
LHTES is the impact of thermal dilatation of PCM, which must be 
taken into account into the conception of experimental and 
operational systems. In both case, geometry differs from the 
idealized version used in most models, with poorly documented 
impact on performances. 


5.4. Analysis: economic impact 


This issue is not treated at all in any of the paper. The use of a 
PCM-based system compared to a standard sensible heat system is 
not discussed. The whole issue of storage in a technico-economical 
context is not found in the reviewed papers. 


5.5. Concluding remarks 


Models are now established for most applications of PCM 
materials. In many cases, simplified models or analytic expressions 
or correlation functions have been developed for practical 
guidance in optimizing design of systems using PCMs. They 
should be used as guidelines when adjusting the implementation 
of a numerical method to be used in the design of systems. 

The validation of numerical predictions should always be 
performed by means of appropriate experimental data. The 
experimental data should in turn be accompanied by a suitable 
uncertainty analysis while the predictions should first address the 
issues of stability, convergence, etc. One should also keep in mind 
that the so-called “governing” equations do not govern but 
describe the physics of the problem. This semantic distinction 
should avoid comments found in a paper into which the authors 
indicate that experiments are not in agreement with the 
predictions instead of the opposite. 

Life cycle analysis, both economical and energetic, should be 
performed on systems to determine into which conditions energy 
storage systems involving PCMs should be used. 

Despite these criticisms, the actual situation which finds 
climate changes, security of the supplies, and fossil fuels depletion 
at the heart of the economic discussions will ensure a place for 
PCMs in the global energy efficiency policies of the world to come. 
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